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The properties of vortex states in a Bose-Einstein condensed cloud of atoms are considered at 
zero temperature. Using both analytical and numerical methods we solve the time-dependent Gross- 
Pitaevskii equation for the case when a cloud of atoms containing a vortex is released from a trap. 
In two dimensions we find the simple result that the time dependence of the cloud radius is given by 
(l + w^t^)^/^ where uj is the trap frequency. We calculate and compare the expansion of the vortex 
I core and the cloud radius for different numbers of particles and interaction strengths, in both two 

and three dimensions, and discuss the circumstances under which vortex states may be observed 
experimentally. 
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I. INTRODUCTION 

(N 

K*" ' Interest in the properties of atomic clouds was greatly stimulated by the experimental discovery of Bose- 

' ' , Einstein condensation in trapped gases of alkali-metal atoms. One of the intriguing features of these novel condensates 

' that remains to be explored experimentally is their behavior under rotation, 
(y-j Vortex states in trapped atomic clouds at zero temperature have been considered theoretically by several authors 

From numerical solutions to the Gross-Pitaevskii equation Dalfovo and Stringari determined the critical 
00 angular velocity, that is the lowest angular velocity for which it is favorable for a vortex to enter the cloud. Lundh et 
On al. Q obtained for reasonably large clouds approximate analytical expressions for the critical angular velocity which 
agree closely with the numerical results. 

For large clouds the presence of a vortex is difficult to detect experimentally, since the size of the vortex core is 
small compared to the size of the cloud. Consequently, the energy of a state in which a vortex is present is nearly 
equal to the energy of the ground state without the vortex. Rather than measuring the total energy of a trapped 
^ ■ cloud it may therefore be advantageous to investigate the density profile of the cloud during a free expansion after the 
Q ' trap potential has been turned off. In this way one may be able to follow how the "hole" in the middle of the cloud 
O ' develops as a function of time, thereby allowing one to distinguish the free expansion of the vortex state from that of 
the ground state without a vortex. In this paper we shall therefore study the free expansion of a cloud containing a 
vortex, which is initially trapped in an anisotropic harmonic-oscillator potential and subsequently released. 

Apart from the difficulties involved in the detection of a vortex state, one also has to consider its experimental 
generation Due to the energy barriers separating rotating and non-rotating states it may be advantageous to 
generate a vortex state by cooling a rotating cloud in its normal state below the transition temperature rather than 
rotating the cloud from rest at temperatures below the Bose-Einstein condensation temperature. 

There are questions regarding the stability of vortices , and the character of the ground state for a fixed angular 
momentum if the angular momentum per particle is not a multiple of ?i However, we shall not address these issues 
in this paper, but will consider the case of expansion of a cloud containing a vortex with one quantum of circulation. 

In the following we shall consider the free expansion of a condensed atomic cloud in the limit when the temperature 
is sufficiently low that the influence of the normal component is negligible. We assume that the system is dilute, in the 
sense that the scattering length is much less than the interparticle distance. In this case one may neglect the depletion 
of the condensate, and the wave function, ^/'(r, t), of the condensed state in an external potential V{r) satisfies the 
time-dependent Gross-Pitaevskii (GP) equation 



ih 



dt 



— V' + Vir) + Uomr)\' 
2m 



(1) 



The effective two-body interaction Uq is given by Uq = ^Trfi^a/m, where a is the scattering length and m is the 
atomic mass. We shall assume the scattering length to be positive, thereby ensuring the stability of the non-rotating 
condensed state for any value of the total number of particles, TV. 
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The plan of the paper is as follows: In Sec. II we consider a two-dimensional vortex state and show that its 
development in time is given approximately by a simple, analytical expression for the cloud radius. We also solve the 
time-dependent Gross-Pitaevskii equation numerically, and compare the resulting density profile to the one obtained 
analytically. In Sec. Ill we discuss the more realistic three-dimensional case by an approximate variational method 
as well as by numerical methods. The expansion of the rotating cloud of atoms is compared to that obtained for a 
non-rotating cloud when the trap potential is turned off. Section IV is a brief conclusion. 



II. VORTICES IN TWO DIMENSIONS 



In the following we consider a two-dimensional geometry, where the vortex is uniform along its axis, which we take 
to be the z-axis. The harmonic-oscillator potential is assumed to be isotropic in the xy-plane and given by 

V = \muj^p\ (2) 

where = + . 

In the stationary case the condensate wave function depends on time only through its phase according to 

^|;{r,t)=',jJ{r)e~'^'''^, (3) 

where ^ is the chemical potential. In the Thomas-Fermi approximation the chemical potential is obtained by inserting 
Eq. (^) into Eq. and neglecting the kinctic-cnergy term. This results in the following expression for the particle 
density: 

IV'P = 7^(A^-^mc.V), (4) 

provided p is less than Pmax — (2/i/mu;^)^/^. For p greater than /Omax; the density is zero. 
The number of particles, v, per unit length along the z-axis is 

z. = 2^ r""dpp|^P = -^. (5) 

In terms of the dimensionless parameter 7, which is defined by 

7 = va, (6) 

the chemical potential is then seen to be given by 

p = Ihw-i^l'^. (7) 

The Thomas- Fermi approximation becomes exact in the limit of large 7. Since the chemical potential is given in 
terms of the total energy per unit length, _E, by the equation p = dE/dv, we conclude that p = 'iE/2v^ and therefore 

E ^ ^vhcj-f^/^. (8) 

Since we shall use Gaussian trial functions later on in discussing the time evolution of the vortex state, it is 
instructive to compare Eqs. and with the result of calculating the ground state energy variationally with a 
trial function of the form 

^j{x,y)^Ae''''/^'' (9) 
where & is a variational parameter. Writing b = aaosc, where a^sc — (h/mujY^'^ is the oscillator length, one finds 

E^.n.i^ + la^ + ^). (10) 

This expression has a minimum for a = aa = {1 + 27)^/^, and therefore the variational estimate of the energy of the 
ground state is 
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E = vhuoy^i-i + l. (11) 

In the limit of large 7 this variational approximation to the energy exceeds the asymptotically exact result (H) by only 
6%. In the opposite limit, for small 7, it is exact to order 7. 

We now consider a vortex state corresponding to one quantum h/m of circulation around the z-axis. The corrections 
to the Thomas-Fermi energy (^) consists of two terms: the vortex energy per unit length, and the kinetic energy 
associated with the rounding of the density profile in the vicinity of the surface of the cloud, which is of the same 
order of magnitude as the vortex energy. For large clouds the vortex energy per unit length is Q 

E. = .|^(0)P^ In (^.^h^^ . ^ ln(3.557^/^), (12) 

where the coherence length ^0 is defined by 



& ; 27V2 



t1 



f7o|V(0)r (13) 



Note that the ratio of the coherence length to the cloud radius in the Thomas-Fermi limit is given by 

^0 fv^ 1 

Pmax 2^ 47I/2 

The kinetic energy associated with the surface is ||] 



(14) 



i?fc = .^ln(2.437^/^). (15) 

Summing up the various contributions, the total energy per unit length for the cloud in the vortex state is therefore 

i?-j.?ic.[^7i/2 + _i_in(13.07)] (16) 

for sufficiently large clouds (7 1). To leading order in 7 the energy of the vortex state is thus the same as that of 
the ground state. We may compare the result (Hq) for the energy to the result of using a variational trial function of 
the form 

V'(p,0)=^pe'^e-'''/2''', (17) 

where is the azimuthal angle, yielding 



E = i^nw\/27-f 4. (18) 



The variational estimate of the vortex energy is obtained by subtracting (11) from (|18[), and gives the value 
i/?ia;3-\/2/47^/^, which is to be compared with Eq. (12). This shows that the variational calculation gives a poor 
estimate of the energy of the vortex for strong coupling. The reason for this is that the trial wave function has only 
one length scale, which determines both the size of the cloud, and the size of the vortex core. 



A. Free expansion 



We shall now use Gaussian trial functions also for treating the time development of the density, when the cloud is 
released from the trap at a certain instant of time, i = 0. We assume the solution to be homologous in the sense that 
the local velocity is proportional to the distance from the axis and therefore employ a trial function of the form 

V'(p, 0, t) = Ape''t'e-p"'^''"e'^P"'\ (19) 

where A — v jixBr^ is the time-dependent normalization constant found from conserving the number of particles 
per unit length. The radial velocity is given by the derivative of the phase with respect to p. The quantities i? and /3 
depend on time. We wish to evaluate the Lagrangian 
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following the method used in Ref. |lO[ . Performing the integration over p, we obtain the Lagrangian as a function of 
the variables P and R and their time derivatives f3 and R: 

From the Lagrange equations for the two independent variables (3 and R we obtain 

mR 



and 



,J « (.3) 



When Eq. for /3 is inserted in (J23|) we arrive at the acceleration equation 



where 



+ (25) 

Note that the contributions to the effective potential energy, U, from the kinetic term due to the zero point motion, 
and the interaction energy both scale in the same way as a function of the size of the cloud. This result is peculiar to 
two dimensions. Equation (^4|) has the solution 

R^{t) = R^{Q)+vlt^, (26) 

where the velocity vq is given by 

Using the result of the minimization in the stationary case with the trial function one finds 

R(Of^aUl + l), (28) 

and thus 

Wo = aoscw(l + |)^/^ (29) 

This implies that the result ( p6| ) may be written in the form 

R{tf = R{Of{l + Loh^). (30) 

This simple result is a consequence of the fact that the effective potential energy varies as R~^, which, as we remarked 
above, is a special feature of two dimensions. 

We note that the root-mean-square (rms) radius prms is V^R, so the final value of the rms velocity is Vrms — V^vq. 
Figure 1 shows how the final value of Wrms varies with 7. For a non-rotating cloud, the corresponding analysis is easily 
carried out with the result tiims = = aosc'jj(l + '^l)'^^^- This is included in Fig. 1. 
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FIG. 1. Final root-mean-square velocity for a two-dimensional cloud in free expansion as a function of the dimensionless 
coupling parameter 7 = i/a. The full line corresponds to the vortex state in the variational treatment and the dashed line to the 
ground state. Numerical solutions to the full time-dependent problem are indicated by crosses (vortex) and squares (ground 
state). The inset shows the time evolution of the rms radius in units Oosc for a vortex (full line) and a ground state cloud 
(dashed line), for 7 = 2.07. 



As we shall see from the numerical solutions presented below, the simple result (^0|) yields an accurate description 
of the expanding vortex state. 



B. Numerical results 

For the numerical study it is convenient to work in terms of scaled quantities, and we introduce the variables 
/ = ij/e^'^{Uo/ nY^"^ , pi — p/aosc and /ii — fi/fiuj. The stationary Gross-Pitaevskii equation may then be written in 
the dimensionless form 



5V 



1^ 

pi dpi 



dpi ^ ^ 



1 

7i 



(31) 



The equation for a non-rotating cloud is similar if the scaled wave function is defined by / = iI)/{Uq/ p)^/"^, and the 
centrifugal term f /2p\ (which comes from the phase e"^ of the wave function in the case of a vortex state) is omitted. 

We integrate this equation using the Runge-Kutta method to find the wave function / as a function of pi corre- 
sponding to a given dimensionless parameter pi. The time-independent wave functions for p = 2.5hu} and p = lOfiio 
are shown in Fig. 2. From normalization, we find that these correspond to 7 = 0.55 and 23.6, respectively. 
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FIG. 2. The radial wave function (solution to the stationary Gross-Pitaevskii equation (^)) for the two-dimensional vortex 
state. The full line is the wave function with chemical potential /ii — 2.5, the dashed line corresponds to fii = 10. Note that 
the units on the y-axis are such that the wave functions are normalized to unity, i. e. they are not the same as in Eq. 



To get a feeling for the orders of magnitude, note that for the realistic case of a system with an axial dimension of 
lOflosc and with Oosc/a = 100, the total number of particles is ~ 10'^7. In the first experiment performed at JILA [Q, 
flosc/i — 700, N — 2000 and the axial dimension of the cloud was of the same order of magnitude as aosc- In the 
loffe-Pritchard trap at MIT ||ll|, Oosc/a = 400, N = 10'^ and the length in the z-direction is SOcosc, giving 7 ~ 30. 
The case 7 = 23.6 is thus quite close to typical experimental conditions. 

When released from the trap, the wave function evolves according to the time-dependent Gross-Pitaevskii equation 
which in dimensionless form is 



dti 2 V opl pi dpi p{ ' 



where ti = ujt. As before, the non-rotating case is obtained by drop pin g the centrifugal term in (p2|). This equation is 
integrated with respect to time using the Crank- Nicholson method |1^ , with the initial value given by the stationary 
wave functions obtained above. During the expansion, we compute the rms radius, Prms, of the system. We find it 
to follow Eq. (H^) to a good degree of accuracy. At large times, pims does indeed grow linearly with time, and the 
final velocities for a few different values of 7 are plotted in Fig. 1. The inset in Fig. 1 shows the evolution of prms 
for 7 = 2.07, corresponding to the chemical potential p, — 3.54?ia; for the vortex state and p = 3.07?iw for the ground 
state. 

The Crank- Nicholson method generates for each time step an explicit table of the real and imaginary parts of the 
wave function. On comparing the wave function at later times with the initial one, we have been able to find that the 
system does indeed undergo a nearly homologous expansion, i.e. the wave function does not change its shape, but 
only flattens out with time. We shall discuss the physical reasons for such behavior, in both two and three dimensions, 
at the end of Section III. 



III. VORTICES IN THREE DIMENSIONS 



A. Simple variational approach 

The clouds studied in experiments do not possess the translational invariance in one direction that was assumed 
in our calculations above. It is therefore of interest to explore the consequences of motion in the third direction on 
the development of a vortex state. Our discussion for the case of two dimensions may easily be generalized to three 
dimensions by employing the trial function 
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(33) 



which describes a cloud undergoing expansion in both the radial and the z-directions. The result is the two coupled 
equations 



R = 



+ 



NUo 



and 



+ 



NUo 



(34) 



(35) 



The inset in Fig. 3 shows the results for the radial and axial expansion (in units ttosc) for the case Na/aosc = 20. 
For large times the velocity becomes constant, just as for the two-dimensional case. In Fig. 3, the aspect ratios of a 
cloud in an isotropic trap before and after expansion are shown as a function of the coupling parameter Na/aosc- The 
aspect ratio, A, is defined as the ratio between the root-mean-square distances, A = Prins/(\/2zi.ms), where the factor 
of \f2 takes into account the fact that there arc two coordinates perpendicular to z. With this normalization, the 
aspect ratio is unity if the rms values of y and z are equal. For the variational wave function we employ, one finds 
Prms = \/2-R and 2;rms = Z, so the aspect ratio is simply R/Z. An aspect ratio differing from unity in an isotropic trap 
would be a clear signal for the presence of a vortex, but from Fig. 3 we find that there are deviations of the aspect 
ratio from unity of more than about 5% only for a number of particles less than lOOoosc/oi- 
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FIG. 3. Aspect ratio for a vortex system in an isotropic iiarmonic trap, before (dasiied line) and after (full line) free expansion. 
The crosses and boxes are computed using the semi-numerical sciicme described in the text. The inset shows the expansion for 
a vortex in the radial (full line) and axial (dashed) directions, for Na/aosc = 20. 



One may ask whether it is possible to improve the signal for the presence of a vortex by changing the external potential. 
We have therefore computed the aspect ratio for a system in an anisotropic external potential V{r) — ^muj'^ [p^ + }? z^) 
as a function of the anisotropy parameter A. An oblate trap corresponds to A > 1, and a prolate one to A < 1. The 
computation was carried out for TVa/oosc = 100, where Oosc = (fi/ {mui\^/'^)Y/'^ is the geometrical mean of the 
oscillator lengths in the spatial directions. Fig. 4 displays the ratio of the aspect ratio of a vortex, Ay, to that of a 
non-rotating system, An, as a function of A, both in the stationary case and after a long expansion time. 
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FIG. 4. Ratio of the aspect ratio for a cloud containing a vortex, Av, to that for a nonrotating cloud, An, in the stationary 
case (full line) and after a long time of expansion (dashed line), plotted as a function of the trap anisotropy A for Na/aosc ~ 100. 
We see that after expansion the difference in the aspect ratios is enhanced for a strong anisotropy, compared to an isotropic 
trap. 

We see that for strong oblate anisotropy, the aspect ratio of a cloud expanding in a vortex state is initially close to 
that for the nonrotating state, but following expansion they are significantly different. For clouds expanding from 
a prolate trap, the ratio of aspect ratios for the vortex and the nonrotating state is enhanced both before and after 
expansion. 

B. Core structure 

In the variational calculations described above, the radial wave function is described in terms of a single length 
scale, R, which determines both the size of the cloud and the size of the vortex core. This is unrealistic when particle 
interactions are important, since the size of the core is then determined by the coherence length, ^, which is a function 
of the density of particles. To study the dynamics of the core we employ a more general trial wave function 



1 



(36) 



with the parameters f3 and Z depending on time. The radial wave function /(p, t), which is normalized to unity, is to 
be determined variationally. 

Minimizing the action obtained from the Lagrangian ( pO| ) with respect to Z, P(t) and f{p,t) yields the coupled set 
of equations 



m 



z = 



m Z 

mz' 

m? Z^ 
2m 



NUn 



(I/O, and 



9p2 pdp 



+ 



NUo 



(I/O 



/ 



I/IV, 



2ttZ 



(37) 
(38) 
(39) 



where (|/|^) = / (Pp\f{p)\'^. These equations may now be solved numerically, at each time step simultaneously taking 
a Crank-Nicholson step for Eq. (|3^) and a Runge-Kutta step for (^). To find the initial conditions, we have to 
consider the static case with an external potential. We choose to consider the expansion of a cloud initially confined 
in an isotropic potential. The relevant equations are obtained from Eqs. (|37[-p9[) by neglecting the time derivatives 
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and, for a potential V{r) = ^muj'^{p'^ + z^), adding to the right-hand side of ( |39| ) the term ^mui'^p'^f and to the 
left-hand side of Eq. ( |38| ) the term uj'^Z. Note that when the particle interaction vanishes, the wave function is that 
of the lowest state with unit angular momentum of a particle in the harmonic oscillator potential. This corresponds 
to putting Z — flosc and f{p) = n^^/^a^^^pe~'' /^"o^o in our variational trial function. Since the Crank-Nicholson 
method preserves the normalization of the wave function and / is to be normalized to unity for any value of the 
coupling NUo, we may produce any stationary wave function of the system by starting out with this one-particle wave 
function and then integrating the static equations of motion discussed above with respect to time, while adiabatically 
increasing the coupling constant up to the desired value p^ . The wave functions we obtain by this method are in 
good semi-quantitative agreement with those obtained by Dalfovo and Stringari [Q from numerical solutions of the 
three-dimensional time-independent Gross-Pitaevskii equation. 

Regarding the time development of the outer radius of the cloud, we find that the results of the semi-numerical 



scheme closely follow those found using the simpler approach, Eqs. (|34| - |35D , especially for weak coupling, as the data 
points in Fig. 3 show. 

We have also considered the time dependence of the core size, which we characterize by the radius, pi, at which the 
particle density first reaches times its maximum value. Results for pi/prms as a function of time are shown in Fig. 

2 

5. If the radial wave function were of the form of the first excited oscillator state, ^ pe~"P , where a is a constant, the 
ratio Pi/pims would be approximately 0.282. This is indeed what we find for Na/aosc ^ 1- For larger values of Na/aosc 
the inner radius is a smaller fraction of Prms initially, but the ratio increases with time. The qualitative behavior 
of the ratio Pi/prms with time for Na/aosc 3> 1 may be understood in terms of two sorts of processes. Initially the 
characteristic time for adjustment of the core size of the vortex. Tad ~ h/nUo, is small compared with the expansion 
time, Tex R/vs, where Vs is the sound velocity of the cloud at t = 0. Under those conditions the core of the vortex 
can adjust essentially instantaneously to the local density, and thus one expects pi to scale as the coherence length. 
Pi cx ^0 oc (R^ /NaY^"^ , or ^ R{R/NaY^^. We take the density n to be N/R{t)"^, since any anisotropy is quite 
unimportant here as long it is of order unity. This behavior ceases at a decoupling time, td, at which Tad = Tox- In 
three dimensions, this means 

R{t) - i?(0)3/4(7Va)i/4 

which gives us 

iotd - (Na/aoscy/^. 

Here we have assumed that R{t) ~ Vgt, which holds when 3> 1. A more careful treatment would result in a 
somewhat larger decoupling time. In two dimensions. Tad — implies that 

R{0) ^ ' 

Since R{t) = R{0)ujt for cut large, the decoupling time is given by 

in two dimensions. At larger times the potential energy will play essentially no role, and the evolution of the cloud will 
be the same as for free particles. To the extent that the cloud is expanding homologously at the decoupling time, we 
expect that pi/prms will remain constant and equal to its value at the decoupling time. The results of the numerical 
calculations exhibited in Fig. 5 are consistent with the assumption of homologous expansion at larger times. 

We conclude that in two dimensions the expansion is homologous to a good approximation during all of the 
expansion, while in three dimensions, the relative size of the core grows until the decoupling time, whereafter it 
approaches a constant value. 

Our results indicate that by allowing a cloud containing a vortex to expand freely, the structure will increase in 
size more rapidly than the size of the cloud. In this way one may thus facilitate optical detection of the structure 
associated with the vortex core. 
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FIG. 5. Evolution in time of the ratio of the vortex core radius to the total system radius for Na/aosc =5 (full line), 10 
(dashed), 20 (dotted) and 30 (dot-dashed), respectively. The straight line shows the value pi/prms =0.282, corresponding to 
the free-particle limit, Na/aosc = 0. 



IV. CONCLUSIONS 

In summary, we have described the time evolution of a freely expanding Bose-Einstein condensed atomic cloud with 
a singly quantized vortex in two and three dimensions using both analytical approximations and direct numerical 
calculation. We have found that simple variational estimates describe very well the time evolution of the radius of 
the cloud. 

In the Umit of weak coupling, the aspect ratio of the cloud with a vortex differs from that of the non-rotating ground 
state both before and after expansion. This effect may be enhanced by using an anisotropic trapping potential. 

In two dimensions, the cloud is seen to undergo a nearly homologous expansion, while in three dimensions, the 
vortex core expands faster than the size of the cloud during the initial stage of expansion, thus making it easier 
to detect experimentally. After a decoupling time td, estimated to be td = lu^^ {Na/aoscY^^ for expansion in three 
dimensions, the size of the core will increase linearly with time. 

We conclude that a vortex is most easily observed experimentally in the weak coupling regime, i. e. for small values 
of the parameter Na/aosc, which is attained either for small numbers of particles or for shallow traps, since the relative 
importance of the vortex to the energy, system size and velocity of expansion is larger in this case. 
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